#### Politics of Pain Script
#### Michael Shepherd
#### 10/26/21

rm(list=ls())
library(readr)
library(haven)
library(foreign)
library(tidyr)
library(dplyr)
library(rdrobust) 
library(rddensity)
library(lmtest)
library(ggplot2)
library(dotwhisker)
library(tidyverse)

###################################################
#       Loading the Data
###################################################

# Border sample of counties within 100 miles of opposing Medicaid expansion border
grdborder <- read.csv("opioid_analyses_data.csv")
# subset of original data for Republican-leaning states 
similarstates <- read.csv("redstates_analyses_data.csv")
# Data on when states expanded Medicaid
medicaidexpansionyear <- read_csv("medicaidexpansionyear.csv")
# Trump tweets about health
tweets <- read.csv("trump_electionyeartweets.csv")
# NYT Opioid Data
nyt <- read.csv("nyt_opioid_mentions.csv")


###################################################
#                 Maps 
###################################################
# Expansion Status Maps for Full and Border Samples

# Prepping data for plots (defining when a case expanded and whether included in border)
medicaidexpansionyear$expanded_15 <- NA
medicaidexpansionyear$expanded_15[medicaidexpansionyear$year_expanded<=2015] <- 1
medicaidexpansionyear$expanded_15[medicaidexpansionyear$year_expanded>=2016] <- 0
medicaidexpansionyear$expanded_15[medicaidexpansionyear$medicaid_expansion==0] <- 0
medicaidexpansionyear$state <- medicaidexpansionyear$state_abv
medicaidexpansionyear$Expansion[medicaidexpansionyear$expanded_15==1] <- "Yes"
medicaidexpansionyear$Expansion[medicaidexpansionyear$expanded_15==0] <- "No"


### Full Sample Map 

plot_usmap(data = medicaidexpansionyear, values = "Expansion")  + 
  scale_fill_manual(name = "Expanded Medicaid", 
                    labels = c("Yes" = "Yes", 
                               "No" = "No"),
                    values = c("Yes" = "#2D2D2D", 
                               "No" = "#D0D0D0"), na.translate=FALSE) + theme(legend.position="right") 


## Border Sample Map 

plot_usmap(include = c("AR", "AZ", "CO", "IA", "ID", "IL", "KS", "KY", "LA", "MD", "ME", "MI",
                       "MN", "MO", "MS", "MT", "ND", "NM", 
                       "NE", "NH", "NV", "OK", "OR", "SD", "TN", "TX", "UT",
                       "VA", "WA", "WI", "WV", "WY"),
           data = medicaidexpansionyear, values = "expanded_15") +  scale_fill_continuous(low="#D0D0D0", high="#2D2D2D",name = "Expanded Medicaid (2015)", label = scales::comma) + 
  theme(legend.position = "none")



########################################
#    Tweets and NYT  plots 
########################################

# Trump Twitter Plot 
p <- ggplot(tweets, aes(x = month_num, y = tweets, color = Topic, group = Topic)) + theme_gray() + 
  geom_line() + geom_point() + ylab("Number of Trump Tweets") + xlab("Month of 2016")

p1 <- p + scale_x_continuous(breaks=seq(0,12,1))
p2 <- p1 + scale_color_manual(values=c("black", "darkgray"))
p2


# NYT Plot 
nyt16 <- nyt[which(nyt$year<=2016),]

x <- ggplot(nyt16, aes(x=year, lty="Number"))  + theme_gray() + 
  geom_line(aes(y=articles))+
  xlab("Year") +
  ylab("NYT Articles")
x2 <- x + scale_x_continuous(breaks=seq(2000, 2016, 4))
x3<- x2 + guides(lty=guide_legend(title=NULL))
x3



##########################################
#    Research Design Assumption Plots
########################################## 

grdborder$Expansion <- as.factor(grdborder$z)

# Pre-trend discontinuity in partisanship? No! 
p <- ggplot(grdborder, aes(x = x, y = dm2008, color = Expansion)) +
  geom_point(size = 0.5, alpha = 0.5) + 
  # Add a line based on a linear model for the people scoring less than 70
  geom_smooth(data = filter(grdborder, x <= 0), method = "lm") +
  # Add a line based on a linear model for the people scoring 70 or more
  geom_smooth(data = filter(grdborder, x > 0), method = "lm") +
  geom_vline(xintercept = 0) +
  labs(x = "Distance to Medicaid Border", y = "Democratic Two-Party Vote (2008)", color = "Expansion")
p1 <- p +   scale_colour_grey()

p1

# Pre-trend discontinuity in opioids? No! 
p <- ggplot(grdborder, aes(x = x, y = opioidrate_10, color = Expansion)) +
  geom_point(size = 0.5, alpha = 0.5) + 
  # Add a line based on a linear model for the people scoring less than 70
  geom_smooth(data = filter(grdborder, x <= 0), method = "lm") +
  # Add a line based on a linear model for the people scoring 70 or more
  geom_smooth(data = filter(grdborder, x > 0), method = "lm") +
  geom_vline(xintercept = 0) +
  labs(x = "Distance to Medicaid Border", y = "Opioid Prescription Rate (2010)", color = "Expansion")
p1 <- p +   scale_colour_grey()

p1


###################################################
#       Effect of ME on Opioid Plot
###################################################

# Impact of Medicaid on Opioids GRD Plot found in the manuscript

grdborder$Expansion2[grdborder$z==1] <- "Yes"
grdborder$Expansion2[grdborder$z==0] <- "No"

grdborder$Expansion2 <- as.factor(grdborder$Expansion2)

p1=ggplot(grdborder, aes(x = x, y = deltaopioids16_14, colour=Expansion2)) 
p2=p1 + stat_smooth(method = lm)
p3=p2+theme_bw()
p4=p3 + scale_color_grey(start=0.4, end=0.2)
p5 <- p4 + ylab("Change in Opioid Rate (2016-2014)") + xlab("Distance to Expansion Border") + labs(color="Expansion") 
p6 <- p5 + geom_vline(xintercept=c(0), linetype="dotted")
p6








###################################################
#            Elections Analyses 
###################################################

# Code used to produce main regression table in the main text 

# Opioid Increase (Regression Table 1, Column 1)
med <- lm(grdborder$shift_16_12 ~  grdborder$increase*grdborder$z + grdborder$x*grdborder$z + 
            I(grdborder$x^2)*grdborder$z  +grdborder$dm2004 +  as.factor(grdborder$stateabbr), weights=grdborder$vap2016)
summary(med)

# heteroskedasticity-consistent standard errors 
med2<- coeftest(med, vcov = vcovHC(med, type = "HC0"))
med2



### Log Model (Regression Table 1, Column 2)

medicaid1_log <- lm(grdborder$shift_16_12 ~  grdborder$log_opioids*grdborder$z + grdborder$x*grdborder$z + grdborder$dm2004 +
                      I(grdborder$x^2)*grdborder$z + as.factor(grdborder$stateabbr), weights=grdborder$vap2016)
summary(medicaid1_log)

# heteroskedasticity-consistent standard errors 
medicaid1_log_b <- coeftest(medicaid1_log, vcov = vcovHC(medicaid1_log, type = "HC0"))
medicaid1_log_b



# Opioid Rate (Regression Table 1, Column 3)

medicaid1 <- lm(grdborder$shift_16_12 ~  grdborder$opioidrate_16*grdborder$z + grdborder$x*grdborder$z + 
                  I(grdborder$x^2)*grdborder$z  +grdborder$dm2004 +  as.factor(grdborder$stateabbr), weights=grdborder$vap2016)
summary(medicaid1)

# heteroskedasticity-consistent standard errors 
medicaid1a <- coeftest(medicaid1, vcov = vcovHC(medicaid1, type = "HC0"))
medicaid1a



################################################################################
#            Special Red state Comparisons 
################################################################################

#### Similar states analyses. Reviewer raised concerns over differences between red and 
# blue states. This analysis relies on a subset of the original border data to just
# Republican-leaning states.

grdborder$Opioid <- grdborder$opioidrate_16
grdborder$MedExpansion <- grdborder$z
similarstates$Opioid <- similarstates$opioidrate_16
similarstates$MedExpansion <- similarstates$z

# full sample
full <- lm(shift_16_12 ~  Opioid + x*MedExpansion + as.factor(stateabbr) + 
             I(x^2)*z  + dm2004, weights=vap2016, data=grdborder)
summary(full)
# results are qualitatively similar in GOP sample
simstate<- lm(shift_16_12 ~  Opioid + x*MedExpansion  + as.factor(stateabbr) + 
                I(x^2)*z  + dm2004, weights=vap2016, data= similarstates)
summary(simstate)


## Coefficient Plot for results 

matrix_coef1 <- as.data.frame(summary(full)$coefficients)  # Extract coefficients in matrix
matrix_coef2 <- as.data.frame(summary(simstate)$coefficients)  # Extract coefficients in matrix

m1 <- NA
m2 <- NA
m1 <- matrix_coef1[c(2,4), ]
m2 <- matrix_coef2[c(2,4), ]

m1$term <- c("Opioid Rate (2016)", "Medicaid Expansion")
m2$term <- c("Opioid Rate (2016)", "Medicaid Expansion")

# Making two standard deviation change for estimates 
2*sd(grdborder$opioidrate_16, na.rm=T)
m1$Estimate[m1$term=="Opioid Rate (2016)"] <- m1$Estimate*85.79484
m1$`Std. Error`[m1$term=="Opioid Rate (2016)"] <- m1$`Std. Error`*85.79484

# no need to multiply expansion as two standard deviations =1 
m2$Estimate[m2$term=="Opioid Rate (2016)"] <- m2$Estimate*85.79484
m2$`Std. Error`[m1$term=="Opioid Rate (2016)"] <- m2$`Std. Error`*85.79484

m1$model<- "Full Sample"
m2$model <- "GOP States"

results_df <- rbind(m1, m2)

results_df$estimate <- results_df$Estimate
results_df$std.error <- results_df$`Std. Error`

dwplot(
  results_df,
  vline = geom_vline(
    xintercept = 0,
    colour = "grey60",
    linetype = 2
  ),
  # plot line at zero _behind_ coefs
  dot_args = list(aes(shape = model)),
  whisker_args = list(aes(linetype = model))
) +
  theme_bw(base_size = 10) + xlab("Predicted Change in Democratic Vote") + ylab("") +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = c(0.66, 0.66),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title.align = .5
  ) +
  scale_colour_grey(
    start = .4,
    end = .6,
    # if start and end same value, use same colour for all models
    name = "Model",
    breaks = c("Full Sample", "GOP States"),
    labels = c("Full Sample", "GOP States")
  ) +
  scale_shape_discrete(
    name = "Model",
    breaks = c("Full Sample", "GOP States"),
    labels = c("Full Sample", "GOP States")
  ) +
  guides(
    shape = guide_legend("Model"), 
    colour = guide_legend("Model")
  )




